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ABSTRACT 



The Hydrographic Sciences Group of the Naval Postgraduate School's 
Department of Oceanography has a goal to establish and continue a research base 
related to the Global Positioning System (GPS). To meet this goal a library of 
programs capable of reducing GPS satellite data to position and time solutions, 
independent of the receiver's navigation solution, is being acquired. 

This thesis describes the modification, documentation, and establishment of 
Texas Instruments 4100 GEOSTAR GPS software written by the Naval Surface 
Weapons Center for their CDC Cyber 865 computer, so that it will run on the IBM 
3033 computer at NPS. Position and time solutions for either a static or dynamic 
receiver can be achieved using this software. 
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I. INTRODUCTION 



The Hydrographic Sciences Group at the Naval Postgraduate School (NPS), 
conducted its first Sea Floor Benchmark Experiment in May 19S5. The experiment 
basically was to deploy a benchmark on the ocean bottom and to establish its position 
with geodetic accuracy. This was to be accomplished using GPS and acoustic range 
data (Kumar, et al., 1984). 

Beyond the goal of the Sea Floor Benchmark Experiment in May 1985 was a 
higher goal of establishing and continuing a GPS research base at NPS. In order to 
reach this goal it has become evident that a library of GPS programs should be 
established at NPS capable of reducing GPS data, independent of outside sources, into 
a format compatible for analysis. 

The Texas Instruments receiver (TI-4100) was selected as one of the GPS 
positioning receivers to acquire data for the experiment. One version of the TI-4100 is 
for commercial users and the other (TI-4100 GEOSTAR) is a tri-agency version for use 
by the Defense Mapping Agency (DMA). National Oceanic and Atmospheric 
Administration (NOAA), and the United States Geological Survey (USGS). TI-4100 
GEOSTAR software has been obtained for use at NPS from NOAA and the Naval 
Surface Weapons Center (NSWC) in Dahlgrcn, Virginia, for compatibility 
modifications with NPS's IBM 3033 mainframe computer. Of this software, three 
programs from NSWC, written chiefly by S. L. Meyerhoff for their CDC Cyber-865 
computer, reformat GPS data and then compute receiver positions and time by either 
an iterative least squares approach or a Kalman filter. 

The main objectives of this thesis are to describe the effort in modifying NSWC's 
software codes and algorithms for compatibility, and establishing them on the 
IBM-3033 mainframe computer at NPS. In addition, discussions of the mathematical 
models, results of test runs and validation of the modified software have been included. 
This software was acquired from NSWC for exclusive use at NTS. Future users at 
NTS, are encouraged to test and verify their solutions as correct before accepting them. 
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II. BACKGROUND 



A. GLOBAL POSITIONING SYSTEM 

GPS is a universal satellite positioning system that is independent of weather, 
time, or geographical position. It provides both position and time in an earth-fixed 
geocentric cartesian coordinate system as well as velocity data. Since it is a space 
based radio navigation system, it can provide accurate information to an unlimited 
number of users globally on or near the earth's surface. Three-dimensional position 
data to within a 16-m spherical error probability (SEP), velocity to within 0.06-0.15 m/s 
and GPS system time to within 25 ns are possible instantaneously to properly equipped 
users anywhere within 500 nmi of the earth (Milliken and Zoller, 1980). 

The GPS is comprised of three major segments, the Space, the Control, and the 

User. 

1. Space Segment 

The GPS when fully operational will consist of 18 active satellites (Space 
Vehicles or SVs) in six planes with near circular orbits at altitudes of 10,900 nmi and 
periods of 12 hours. Each of the six orbital planes is inclined 55 degrees to the equator 
and is to contain three SVs equally spaced. The configuration of the SVs in the orbital 
planes is such that a minimum of four SVs will be available to the user at all times, 
anywhere on or near the earth's surface. Each individual SV transmits its broadcast 
ephemeris and clock correction, plus an almanac which contains orbital parameters and 
clock correction estimates for all other SVs in the system. 

Each SV transmits two radio frequencies simultaneously, LI and L2 
respectively, containing navigation information in the navigation data message. The 
LI frequency is centered on 1575.42 MHz and the L2 frequency is centered on 1227.60 
MHz. These signals are in a frequency range that provides good all weather operation 
and allows the user to determine ionospheric propagation delay. The LI signal is 
modulated by a 10.23 MHz clock rate Precision (P) code and a 1.023 MHz 
Course/Acquisition (C/A) code (Spilker, 1980). The P-code is a random sequence 
binary code that provides the user the capability for high-precision positioning and is 
highly resistant to electronic and multipath interference. The C/A-code gives all users 
rapid acquisition capability and acts as an aid to gain access to the P-code. The 
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navigation data message allows the user's receiver to calculate position, velocity and 
GPS time by providing broadcast ephemeris, clock corrections and almanacs for all 
SVs. The L2 signal is modulated by the P-code only. 

2. Control Segment 

The control segment consists of a master control station, a number of monitor 
stations and ground antennae located around the world. The USAF master control 
station is at Colorado Springs, Co., and monitor stations are located on Kwajalein 
Atoll, Ascension Island, Hawaii and Diego Garcia. DMA also maintains three stations 
in the United Kingdom, Australia and Argentina, and stations are proposed for 
Ecuador and Bahrain. The monitor stations acquire ranging data for each SV in a 
passive mode. These data are then relayed to the master control station where they 
are processed to determine the SV's position and signal data accuracy (Wooden, 19S5). 
The master control station takes the relayed data, computes errors, generates a new 
navigation data message and then uploads this message to each SV via the ground 
antennae. The ground antennae transmit as well as receive satellite control 
information. 

3. User Segment 

This segment consists of stationary, low, medium and high dynamic receivers 
designed with specific requirements for individual users. The user receivers are 
designed to receive and process SV data from four SVs either simultaneously or 
sequentially. The TI-4100 receiver determines the range from the SV by measuring the 
time delay of the SV's specific epoch and the user generated codes. Then making 
separate ranging and carrier frequency Doppler shift measurements to four SVs, the 
receiver calculates the user's position, clock bias error, clock drift and velocity. The 
position coordinates are in the World Geodetic System (WGS) 1972, an earth-fixed 
geocentric cartesian coordinate system. 

B. TEXAS INSTRUMENTS 4100 NAVSTAR RECEIVER (TI-4100) 

1. Characteristics and Description 

The TI-4100 receiver is portable for field operations (it is about the size of a 
small suitcase) but can be mounted in a standard 19-inch rack if desired. A small 
lightweight hand-held Control Display Unit, contains a keyboard and a display 
window. The receiver is modular in design and allows for easy removal and 
replacement of circuit boards. 
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The Tl-4100 is a single channel dual-frequency, multiplexing NAVSTAR 
receiver and navigation processor that tracks up to four SVs (Texas Instruments, 
1983). The dual frequency receiver allows for tracking LI and L2 frequencies for each 
of the four SVs. It also has an optional dual cassette recorder for recording GPS/SV 
data. If GPS/SV data is to be acquired in a degraded mode (i.e. less than four SVs) 
there is the capability of adding an external atomic clock time standard input to the 
Tl-4100 receiver. The Tl-4100 receiver, Control Display Unit, antenna and cassette 
recorder are shown in Figure 2.1. 

2. Receiver Measurements and Postion Calculations 

The basic measurements made by the Tl-4100 receiver are the arrival times of 
the SV-generated signals compared to the receiver clock, which contains a stable 
oscillator and a counter (Texas Instruments, 1983). Four pseudorange measurements 
plus the navigation data message provide enough information to calculate a 
three-dimensional user position solution plus the user clock correction to true GPS 
time. The solutions are obtained using four observation equations and four unknowns. 
Doppler shifts due to the relative velocity between the SV and the user are measured. 
Measurements are also made of the phase differences for the LI and L2 signals 
between the SV and the user. 

3. Position Accuracy 

Error statistics for the Tl-4100 receiver are given in Table I. It should be 
noted for this table that Operating Modes using time aiding, have a positive growth 
rate error of 2.0 to 10.0 m/h (1-sigma) for frequency standards that allow navigation 
calculations to estimate precise frequency (Texas Instruments, 1983). The positioning 
accuracies in Table I were based on the Tl-4100 receiver observing pscudoranges, 
Doppler shift data, phase differences for LI and L2, then computing position, time and 
velocity solutions. The pseudoranges from the Tl-4100 receivers are the major 
observables considered in this thesis. 

C. PSEUDORANGE DATA AND TI1E1R CHARACTERISTICS 

The basic operation of the GPS depends on the user determining pseudorangc 
and range rate to a number of GPS SVs which have precisely known ephemerides 
(Texas Instruments, 1983). Pscudorange is the apparent range, which includes the 
clock errors between the SV and the user's receiver. The transit time of the navigation 
signal between each SV and the user is measured by the receiver and then scaled by the 
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Figure 2.1 Tl-4100 NAVSTAR System 

(From lexas Instruments, 1 1-4 100 Owner's "Manual, 1983). 
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TABLE I 

TI-4100 POSITION ACCURACY 



(From Texas Instruments, TI-4100 Owner's Manual, 1983) 



Operating Mode 


Error 

Stati sties 


1-Sigma 

Position 

p ("0 


i 

Error 

C/A 


1-Sigma 
Veloci ty 
Error 
(m/s) 


Four satellites 












Stati onary 

0. 5 g acceleration 

2 g acceleration 

4 q acceleration 

Di f ferential 

Stationary differential 


3-D 

3-D 

3-D 

3-D 

3-D 

3-D 


13 

14 
14 
30 

4 

3 




46 

47 
47 
55 
12 

8 


N/A 
0. 15 
2. 00 
5. 00 
0. 15 
N/A 


Three satellites 












Degraded (altitude aiding) 


Hori zontal 


16 




24 


0. 15 


Differential degraded 
(altitude aiding) 


Hori zontal 


5 




6 


0. 15 


Degraded (external cesium 
frequency standard time 
aidi ng) 


3-D 


16 




47 


0. 15 


Differential Degraded 

(external cesium frequency 
standard time aiding) 


3-D 


5 




12 


0. 15 


Two satellites 












Degraded (altitude and 
time aiding) 


Hori zontal 


16 




24 


0. 15 


Differential degraded 

(altitude and time aiding) 


Horizontal 


5 




6 


0. 15 



speed of light to compute a raw pseudorange. Each raw pseudorangc must be 
corrected for tropospheric and ionospheric propagation delays, SV clock olTset, user 
clock bias, relativity corrections, and earth rotation corrections. The pseudorangc 
calculation is given by Equation (2.1). 

R i* = R i + cAt Ai + c(At u -At si ) (2.1) 

where i = 1 to 4 (pertaining to four observed SV's) 
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Rj* = pseudorange to the SV 
Rj = true range 
c = the speed oflight 

At<- = SVj clock offset from GPS system time 
At u = user clock offset from GPS system time 
AtAj = propagation delays and other errors 

The tropospheric delay is estimated using simple refraction models with standard 
atmospheric values. Ionospheric delay is calculated by taking advantage of the dual 
frequency measurements of the LI and L2 signals. The SV clock error (difference 
between SV clock and GPS system time) is modeled by the Control Segment and these 
correction coefficients are transmitted in the navigation data message from the SV to 
the user. User clock bias is due to the fact that the receiver clock is not perfectly 
synchronized to the SV clock. Relativity corrections are due to the eccentricity of the 
SV's orbit. Earth rotation corrections are due to the motion of the earth during the 
signal's propagation from SV to the user. Once four raw pseudoranges have been 
corrected for systematic errors, they are used to calculate the receiver's 
three-dimensional position and time. 

The precision of position calculated by pseudoranges is dependent on the 
geometry of the SV's with respect to the user. This geometry is continually changing, 
even if the user's receiver is not in motion relative to the earth. If the intersection of 
the pseudoranges is orthogonal, errors in the pscudorange measurement arc minimized 
whereas if the angle of intersection is small, errors are large because of the geometric 
relationship. 

Geometric Dilution of Precision (GDOP) is a composite "indicator" expressing 
the influence of the SV's geometry on the accuracy of the user's clock oflset and 
position determinations (Milliken and Zoller, 19S0). Other terms to describe dilution 
of precision in position and time are: Position Dilution of Precision (PDOP) which 
yields dilution of precision in three-dimensions; Horizontal Dilution of Precision 
(HDOP) which gives dilution of precision in the two horizontal dimensions; Vertical 
Dilution of Precision (VDOP) which yields the dilution of precision in the vertical 
dimension; and finally Time Dilution of Precision (TDOP) is the dilution of precision in 
the user's time bias (Milliken and Zoller, 1980). GDOP, which is dimensionless, is 
calculated by Equation (2.2), (Milliken and Zoller, 1980), 
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( 2 . 2 ) 



GDOP = [(PDOP) 2 + (TDOP) 2 ] 1 / 2 

The best geometry for a four SV configuration is the combination which results in the 
lowest GDOP. 

Once pseudorange, Doppler shift and phase data have been acquired, and the 
receiver's field solutions calculated, post-processing of these data may yield improved 
position and time solutions. To make TI-4100 data accessible to other users with 
different computer processing systems, it has become necessary to format SV data in a 
manner that can be exchanged among independent users. 

D. STANDARDIZED EXCHANGE FORMAT FOR NAVSTAR GPS 

GEODETIC DATA 

This Standardized Exchange Format (SEF) was proposed by V. Dan Scott and J. 
Gary Peters (1983), and was adopted to make possible the universal exchange of 
NAVSTAR GPS data. To do this, serious compromises were made in regard to 
optimum magnetic tape utilization, computer execution efficiency, and computer 
program size (Scott and Peters, 1983). These compromises were required to allow for 
maximum tape compatibilty between different computer systems. 

1. The Basic Format 

The building block for each standardized data exchange tape is the record, 
which is 80 characters. Under this format, an unlabeled 9-track, ASCII-coded tape 
with an 80-character blocking factor is used. The header (or the first) file on a tape 
consists of comment records and is not considered part of the data set but is used to 
help the user to identify the contents of the tape. When the last data file is written on 
the tape a special end of information file is written (Scott and Peters, 1983). 

On the tape, records fall into two categories: Control Information or Data 
Information records. Control Information records are used to identify which 
FORTRAN Format statement to use when reading data records. Data Information 
records consist of one of eight data items that are uniquely identified by a four digit 
data item number. The eight data items under the SEF are listed below: 

(1) Constant equipment and software specific data 

(2) Campaign and station specific data 

(3) Calibration data 

(4) Time tag update data 

(5) Space Vehicle specific data 
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(6) Weather data 

(7) Real time solution data 

(8) Measurement data 

The goal of this format is to be independent of hardware and software design. In 
addition to data formats, the reference Coordinate System for the data must be defined. 
2. Coordinate System 

The WGS 72 an earth-fixed geocentric cartesian coordinate system is defined 
as follows (Seppelin, 1974): 

X-Axis = Intersection of the WGS 72 reference meridian plane 
and the plane of the Equator. 

Y-Axis = Measured in the plane of the Equator 90 degrees 
east of the X-Axis. 

Z-Axis = Parallel to the direction of the conventional 

international origin (CIO) for polar motion. 

The reference ellipsoid for WGS 72 is defined by the parameters semi-major axis, a = 
6,378,135 m, and flattening, f = 1/298.26 (Seppelin, 1974). 
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III. DATA PROCESSING 



The pre-processing of GPS data is carried out to obtain a format suitable for 
further manipulation. This is necessary when geodetic positioning of the receiver 
(either static or dynamic) is desired, through post-processing the SV data. Figure 3.1 
portrays the flowchart for processing the TI-4100 data. 

A. PRE-PROCESSING 

1. Manipulation of Cassette Recorded Data 

The first step in processing GPS data recorded on cassettes, is to dump the 
raw data from the cassettes to either a 9-track magnetic tape or disk. The Applied 
Research Laboratories (ARL) at Austin, Texas have developed software to accomplish 
this fundamental task. The software at ARL can dump and process the raw data 
recorded on cassettes, and format the GPS data onto a 9-track magnetic tape in the 
SEF (Chapter II. D.). This tape can then be transferred to disk for convenience in 
further processing with a mainframe computer. 

2. Alignment of Common Time-Tagged Data 

The alignment of common time-tagged data items is accomplished through use 
of the CON9TR program. 

a. Program Description 

Program CON9TR, the first of three programs acquired from NSWC 
(Chapter I.), reads the SEF data and formats it into another NSWC specified file. This 
file then becomes the input for the two remaining programs, SOLTOM and KALMN2 
and has all common time-tagged data on the same record. The only calculations 
performed at this stage are variances on the pseudorange data and the delta range 
(Doppler) data which are measures of the signal to noise ratio. 

The program also looks at the "quality vectors" (Texas Instruments, 
1982) of the data on the 9-track SEF tape and returns an integer value representative 
of the data quality at that time line. Four integer values are possible: 

0 = good data 

1 = good data but this is the first Doppler count 

14 = bad data at that point but there was not a break 

in the Doppler data 

15 = bad data 
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Figure 3.1 Flowchart For Processing TI-4100 GFOSTAR GPS data. 
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The CON9TR program requires a user input file (Figure 3.1) in addition 
to the GPS data in the Standardized Exchange Format during program execution. 
CON9TR consists of the main program and 22 subroutines. Appendix A gives a brief 
description of each subroutine's function and nesting. 

b. Input/Output 

The user input file (Figure 3.1) for the CON9TR program consists of the 

following: 

(1) The number of trackers, usually four but may be less than four 

(2) GPS time of the first data point to be put on the new file 

(3) Time to end data on the new file, if the end of the file 

is found before the time to end data, then data will end 
at the end of file 

(4) Range variance bias factor [km ] 

(5) Doppler variance bias factor [km ] 

(6) X-coordinate of the initial receiver position that 

will be put on the new files 

(7) Y-coordinate of the initial receiver position that 

will be put on the new files 

(8) Z-coordinate of the initial receiver postion that 

will br put on the new files 

Appendix B gives the user input guide to CON9TR program variables and format. 

The first output file from CON9TR (Figure 3.1) writes all common time 
satellite tracking data into the same record block, and Appendix C gives the output 
format. The second output file writes out messages to aid the user in data quality 
analysis, flagging one of four quality vector conditions that exist in the data. 

c. Compution of Variances 

The variances for the pseudorange and Doppler data in the CON9TR 
program are computed by Equations (3.1) and (3.2), respectively (Texas Instruments, 
1982). 

a2 code = l(B L *N)/(10( ( V 10 ))]* (3.1) 

[0.5 + Bp O /(10< C i/ 10 >)] * F m + K 2 

where the asterick is used to signify multiplication 
Bj = code loop noise bandwidth Hz 
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Bpo = code loop predetection bandwidth Hz 
C; = Lj carrier power to noise spectral density ratio 
dB-Hz 

N = tracking mode, 1, 2, 3, 4 based on SV mode 

9 9 

Kj = (chips) - to [m - ] conversion: 

(29.305) 2 for P-code 
(293. 05) - for C/A-code 
K-> = variance bias factor [m 2 ] 

ff2 pLL = K B L * X )/( lo(C i 10) )1 * 0 + [0.5 * (3.2) 

b pd (io( c i 10 ))]} * k 3 + k 4 

where, 

PLL = carrier phase-locked loop 
B p = carrier loop bandwidth Hz 
Bppj = carrier predetection bandwidth Hz 

9 

K 3 = [rad - ] to [m - ] conversion 
= (c/27tf) 2 
c = speed of light 
f = L band carrier frequency 
Lj = 1.57524 x 10 9 Hz 
L 2 = 1.2276 x 10 9 Hz 

“ 9 

K 4 = variance bias factor [m - ] 

d. Program Modifications at NPS 

The original code from NSWC was written in FORTRAN IV. This 
code was modified at NPS to be compatible with the FORTRAN 77 compilers 
(FORTVS and WF77) on the IBM-3033 mainframe computer. Changes to the original 
code included the following, 

• Changes in the formatted write statements 

• Changes in variable alignments in the common statements 

(i.e. making sure double precision variables are aligned 
ahead of single precision variables) 

• Changes in dimension size of the data arrays 

• Changes in defining variables 
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• Changes in the placement of data intialization 

statements 

• Addition of the Block Data subroutine to contain the 

moved data intialization statements and accompanying 
common statements. 

Executive files, written at NPS to handle the input and output files for 
the CON9TR program, included the auto-double precision commands for all real 
variables using the FORTVS compiler. This modification automatically doubles the 
declared precision in the code and was required due to the use of 60 bit words by the 
CDC Cyber-865 computer and the 32 bit words by the IBM-3033 mainframe computer 
at NPS. 

B. DETERMINATION OF PSEUDORANGE POSITION SOLUTION 

User receiver position solutions whether fixed (static) or moving (dynamic) may 
be calculated by post-processing the TI-4100 GEOSTAR GPS data. 

1. Static Positioning 

The second program, SOLTOM from NSWC (Chapter I.) computes a position 
solution for a static receiver by using an iterative batch least squares calculation. 

a. Program Description 

Program SOLTOM computes receiver position and time, solving for 
initial corrections to estimates of parameters X, Y, Z, time biases and scaling factor for 
a tropospheric correction. Types of solutions which can be selected through the user 
input file are described in Chapter lII.B.l.e.. The program consists of the main 
program and 52 subroutines. Appendix D gives a brief description of each subroutine's 
function and nesting. 

b. Math Models 

The math models for the batch least squares solution in program 
SOLTOM solve for corrections to the initial estimates of the parameters desired. The 
observation equation is given by Equation (3.3): 

R°i = [(X sy - X p ) 2 + (Y SV - Y p ) 2 + (Z $v - Z p ) 2 J 1/2 + CB + CD + CA + TR (3.3) 



where, 



R°j = i th 



range observed from SV of receiver 
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X sv = earth-fixed geocentric X-coordinate of SV 

Y sv = earth-fixed geocentric Y-coordinate of SV 

Z $v = earth-fixed geocentric Z-coordinate of SV 
Xp = earth-fixed geocentric X-coordinate of receiver 
Yp = earth-fixed geocentric Y-coordinate of receiver 
Zp = earth-fixed geocentric Z-coordinate of receiver 

CB = Clock bias term, SV or receiver depending on user 
input selection 

CD = Clock drift term, SV or receiver depending on user 
input selection 

CA = Clock aging term (frequency drift), SV or receiver 
depending on user input selection 

TR = Tropospheric refraction correction. 

The least squares solution of Equation (3.3) in matrix form, can be written as: 

BP + E = 0 (3.4) 



where, 

B = Matrix containing the partial derivatives of the observation 
equation with respect to the parameters being solved for 

P = Vector containing the differences of the parameter adjusted 
(Pa) and the parameter estimate (Pe), (Pa - Pe) 

E = Vector containing the partial derivatives of the observation 
equations multiplied by the weight matrix, multiplied by the 
data residuals 

The matrix B is obtained from Equation (3.5), 

B = A T wA (3.5) 
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where, 

A = Partial derivatives of data equation with respect to a 
parameter j 

= Transpose of the A-matrix 

J 

w = Weight matrix with a ^ on the diagonal, zeros 

off diagonal, where k = 1 to number of data points, this 
is the variance of noise level on a "k" range measurement 

or writing Equation 3.5 in summation form, 

B i,i - IV i aR \ / ap i * * V o.6) 



where, 



n = number of data points used in the least squares solution 

R^o = k^ data equation 

Pj = i 1 ^ parameter being solved for 

Pj = parameter being solved for 

0^2 = k^ variance value 

The A-matrix or design matrix is defined by Equation (3.7): 




where, 

r°. = jth c j ata p 0 j nt (33) 

Pj = j n parameter being solved for 

The inverse of B is the variance-covariance matrix, COV, with the variances <7~ of the 
parameter corrections AP that arc being solved for, on the diagonal given by Equation 
(3.8), 
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B' 1 = COV 



( 3-S > 



where, COV is a 6 x 6 matrix if six parameters are being solved, or “ x ” if the optional 
7 Ln parameter (CA) is to be determined. The off diagonal elements are the 
cross-correlations of the paired parameters for that time line. The vector E in 
summation form is shown in Equation (3.9). 



E; = 



IVi^k 



°c 



p i> 



R-, 



es 



(3.9) 



where. 

R-m = k data equation 

n = number of data points used in least squares 
calculation 

P- = ri" parameter being solved for 
<7-^~ = k l “ weighted variance value 
R rs$ = range residuals 
or writing Equation (3.9) in matrix form. 



E = R. es wA (3.10) 

where. 

AP = B^E (3.11) 



and 

AP = corrections to the initial estimates of unknown parameters being solved for (X 
Y, Z. clock biases, and tropospheric refraction correction i. 

The following Equations (3.12) to (3.35) describe the various corrections. Time Epoch, 
clock. Ionospheric refraction. Tropospheric refraction, earth rotation, and relativistic 
effects which are required to correct pscudorange and delta range (Doppler) data 
(MeyerhofT. 19S5). 






Time Epoch Correction 




(3.12) 



where, 

T $y — time from epoch 

tt = the GPS time-tag for range and delta range (Doppler) 



T qc = time of reference for clock corrections 

To correct the T $y for the end of the week cross over, correct T sy as follows: 

If T sv < 302400 [s] then, T sy = T $y + 604800 [s] 

If T sy > 302400 [s] then, T sy = T $y - 604800 [s] 



where, 

C c ](tt) = clock correction for pseudorange data 
at GPS time-tag in [km] 

AC c j(tt) = clock correction for delta range 

(Doppler) data at GPS time-tag in [km] 

A Q = clock correction bias term [s] 

Aj = clock correction drift term [s/s] 

Aj = clock correction aging term 



data 



Clock Correction 



C d («) = IA 0 + (Ai * T sy ) + (A 2 * T sy - T sy )] * c 



(3.13) 



AC d (tt) = [A 0 + (Aj * {(2 * T sy ) - td})] * td - c 



(3.14) 



]s/(s) 2 ] 



T sy = time from epoch 
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c = speed of light 299792.458 [km/s] 

td = the Doppler count interval for delta range 
(Doppler) data 

Ionospheric Refraction Correction 

C iQ (tt) = CPj(tt) - CP 2 (tt)/(Q 1 /Q 2 ) 2 - 1.0 (3.15) 



DCjo(tt) = [(POP T !(tt)/Q,) - (DOPj ? (tt)/Q ? )] * c (3.16) 

[(Qp'Q 2 ) 2 -L0]V os 

where, 

Cj 0 (tt) = pseudorange ionospheric refraction 
delay in [km] 

tt = GPS time-tag for pseudorange and delta range 
(Doppler) 

CPj = LI pseudorange in [km] 

CP 2 = L2 pseudorange in [km] 

DOPpi = LI Doppler count 

DOPl 2 = L2 Doppler count 

DCj 0 = delta range (Doppler) frequency shift due 
to ionosphere 

Qj = LI frequency multiplier (154.0) 

Q 2 = L2 frequency multiplier (120.0) 

V QS = The nominal satellite oscillator frequency, 10.23 MHz 

In Equation (3.15) the correction, Cj 0 (tt), is used to correct the pseudorange that is 
time- tagged, tt. In Equation (3.16) the correction, DCj 0 (tt), is applied to the delta 
range (Doppler) data to correct it for time-tag, tt. 
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Tropospheric Refraction Correction 

The Tropospheric refraction computation can be divided into two parts, 
weather dependent and receiver geometry calculations (Meyerhoff, 1985). 

Weather Dependent Correction 

First, convert temperature to degrees Kelvin and convert humidity to a 

fraction: 

T k = T c + 273.0 (3.17) 



RH = H d / 100.0 (3.18) 

where, 

T £ = air temperature in degrees centigrade 
H d = the relative humidity in percent. 

Second, compute the surface partial pressure, E Q , from temperature and relative 
humidity. 

E 0 = RH * 35.65 * 10 exp[7.617-(2285.0/T k )] (3.19) 

Third, the components of the zenith value of tropospheric delay, Z dr y and Z wet are 
calculated from, E Q , and atmospheric pressure. 

Z dry = 2.276 * Pb * 0.01 (3.20) 



Z wet = ( 47 ° * E 0 exp(1.23)/T k ) + (1.705 * 10 6 (3.21) 

* A h " E 0 cxp(1.46/T k 3 ) 

where, 

exp = exponent 

Pb = the atmospheric pressure in [mb] 

= temperature lapse rate set at 0.006 °C/m 
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exp — exponent It should be noted that Z ^ . and Z wet need only be recomputed, 
when the weather data changes. 

Receiver Geometry Dependent Correction 

First, to find the elevation angle, EV(t), of the satellite at time, t, take 
the inverse sine of the dot product of the unit vector from the receiver to the satellite 
and the unit position vector of the receiver. 

EV(t) = sin' 1 R$T 1? (3.22) 

where, 

RS = vector from the receiver to the satellite 
MG rs = vector's magnitude 
R = position vector of the receiver 
MG r = magnitude of R 



MG rs = [(X(t) - X r ) 2 + (Y(t) - Y f ) 2 + (Z(t) - Zj.) 2 ] 1 / 2 (3.23) 

MG f = [X r 2 + Y r 2 + Zj. 2 ] 1 / 2 (3.24) 

where, 

X f , Y f , Z r = the estimate of the receiver's 
position in an earth-fixed 
geocentric coordinate system 

rewriting Equation (3.23), EV(t) can be written as, 

SEV(t) = K(X(t) - X r ) * X r ) (3.25) 

+ {(Y(t) - Y r ) * Y r )} 

+ {(Z(t) - Zj.) * Z r }]/(MG rs * MG r) 

where, 

X(t), Y(t), Z(t) = the position coordinates in a earth- 
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fixed geocentric coordinate system 



or 



EV(t) = sin' 1 [SEV(t)] (3.26) 

Second, the multiplication factors for dry and wet as a function of elevation are 
calculated from the elevation angle, EV,at GPS time-tag (tt). 

F dry (t) = 1 / [sin EV(t)J + {0.00 143/K tan EV(t)) + 0.0445]} (3.27) 



F wet (t) = 1 / [sin EV(t)] + {0.00035/[(tan EV(t))+ 0.17]} (3.28) 

Third, F dr y(t - td) and F wet (t - td) are calculated in a similar manner as Equations 
(3.27 and 3.28). 

Total Tropospheric Refraction Correction 

TR(tt) = [Z dry * F dry (t)] + [Z wet * F wet (t)] (3.29) 



DTR(tt) = {Z dry * [F dry (t) - F dry (t - td)]} + {Z wet * [F wet (t) - (3.30) 

W* - td )» 

where, 

TR(tt) and DTR(tt) are range corrections are in [km]. 

Earth Rotation Correction 

An earth rotation correction is required due to the rotation of the earth 
during the signal propagation. 

C er = l(Y(t) * X r ) - (X(t) * Y r )j * RRE/c (3.31) 
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DC er = [({Y(t) - Y(t - td)} * X r ) - ((X(t)-X(t - td» * Y r )J * RRE/c 



(3.32) 



where, 

X(t), Y(t), Z(t) = position coordinates in an earth- 



fixed geocentric cartesian 
coordinate system at pseudorange 
transmit time "t" used with data 
time tagged "tt" 



X p Y f = estimate of the X and Y coordinate 



in an earth-fixed geocentric 
cartesian coordinate system 



RRE = Rotational rate of the Earth corrected for motion 
of the equinox 
(0.00007292115855 [rads/s] 

td = Doppler count interval for delta range 
(Doppler) data. 

Relativistic Correction 

The correction for the relativistic effects is due to the eccentricity of the 
satellite's orbit. 



Rel(tt) = (2 * [GMj'^/c) * e * A * sin [Ek(t>] 



(3.33) 



at the beginning of the Doppler count, 



DRel(tt) = (2 * [GM]^ 2 /c) * e * A * sin [Ek(t - td)] 



(3.34) 



W’hich yields, 



DRC(tt) = RC(t) - RC(t -td) 



(3.35) 



where, 

tt = GPS time-tag for both pseudorange and delta range 



(Doppler) data 
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td = Doppler count interval for delta range (Doppler) 
data 

t = GPS time-tag at transmission time for data time- 
tagged tt 

e = eccentricity of the satellite orbit 

A = the square root of the semi-major axis of the 
satellite's orbit 

Ek(t) = eccentric anomaly for the satellite at pseudorange 
transmission time t 

Ek(t - td)= eccentric anomaly for the satellite at the 
beginning of the Doppler count interval 

GM = WGS 72 value of the Earth's gravitational constant including 
atmosphere 39S600.8 km^/s* 

RC = computed range 

these are the defined constants and variables. 

c. Input/Output 

There are two input files for the Soltom program (Figure 3.1). The first 
file is the user input file, where the user sets which types of solutions will be calculated 
and which data solutions will be printed out. The user input guide for this file is 
described in Appendix E. The second input file for SOLTOM is an output file from 
the CON9TR program containing records that are aligned with common time-tagged 
GPS data items. 

There are also two output files for SOLTOM (Figure 3.1). The first file 
is the updated receiver position, time and other print options such as solution sigmas 
and clock bias terms, (as determined by the user input file). The second output file is a 
data corrections residual file, which can be used in plotting routines and for input to 
other software that handle multiple station positioning (Meyerhoff, 1986). 

d. Program Modifications at NPS 

The source code from NSWC for the SOLTOM program was written in FORTRAN 
IV which has been modified to run at NPS on the IBM-3033 mainframe computer. 
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Similar changes to the source code in the CON9TR program were also made to the 
SOLTOM program including the executive files (Chapter I1I.A.2.). The original code 
was written to be able to handle the precise ephemeris data as well as the onboard 
broadcast ephemeris data in the navigation data message from the SV. The test data 
for SOLTOM was generated at XSWC using the Broadcast ephemeris only, and as 
such, the modified program version at XPS has the Precise Ephemeris subroutines and 
corresponding, call, statements commented out. 

e. Solution Techniques Through Least Squares 

There are ten types of solutions possible with the SOLTOM program 
which are dependent on the user input file (Meyerhoff, 1985). Appendix E gives the 
details for the type of solution desired. Of these ten types, nine are batch least squares 
solutions and one solution type five is a Xewton-Raphson solution (a geometrical 
solution, not a least squares fit). The type-five solution can be used for either static or 
dynamic positioning performing an independent solution at each time line. The nine 
batch least squares solution types fall into three data classes, Range (pseudorange), 
Range Difference (delta range Doppler), and Bias Range. 

The reader is WARNED, that the Bias Range (solution type ten) has still 
not been tested by NSWC, and thus it should not be used at XPS. Solution types one, 
four, six, and nine are Range (pseudorange) solutions. Solution types two and seven 
are Range Difference (Doppler) solutions. Finally, solution types three and eight are 
combination solutions of the Range (pseudorange) and Range Difference (Doppler) 
data. 

The SOLTOM program has the capability to run more than one 
iteration of batch least squares, if needed (e.g. for bad initial receiver position or time 
bias problems). The user has control via the user input file to set how often the least 
squares solution will be calculated and printed out. When a solution is to be 
computed, the program takes the B-matrix containing the partial derivatives of data 
equations with respect to the parameters being solved for. and multiplies by the 
variances of the range measurements (Equations 3.5 and 3.6), inverts the B-matrix, 
which is the variance-covariance matrix with variances of the corrections AP on the 
diagonal, and computes the corrections to the initial estimates of the parameters being 
sought, AP (Equation 3.11). Then Equation (3.36) gives the updated station position 
coordinates, clock bias, drift, aging, and Tropospheric refraction correction: 

X = X + AP X (3.36) 
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Y = Y + APy 



Z = Z + AP Z 

A o — A o + A ^Ao 

A l l = A l l + AI> Ajt 
A 2 t 2 = A 2 t 2 + AP A2t 

TR = TR *b APjj^ 
where, 

X, Y, Z = updated earth-fixed geocentric cartesian 
coordinates through batch least squares 

A q = clock correction bias term in [s] 

Aj = clock correction drift term in [s/s] 

A 2 = clock correction aging term in 
[s/(s) 2 ] 

(optional 7th parameter to solve for) 

TR = tropospheric refraction correction 

The new solution takes all previous solutions and averages up to the current time line, 
before producing the iterated solution. 

2. Dynamic Positioning 

The third program, KALMN2 from NSWC computes a position solution for a 
moving (dynamic) receiver, using the UDU' factorization of an eight-state extended 
Kalman filter (Chapter I.). 

a. Program Description 

Program KALMN2 uses the first four-states for receiver position (X, Y, 
Z) and receiver clock bias, and the other four-states for receiver velocity along three 
axes (X, Y, Z) and receiver clock drift. By using the Doppler phase (change in range) 
measurements to smooth the pseudoranges, signal multipath and measurement noise 
effects can be suppressed significantly to produce a viable solution (Meyerhof!' and 
Evans, 1986). 
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Note: The four-states of the velocity and clock drift arc still under 
testing at NSWC. 

The KALMN2 program consists of the main program and 43 
subroutines. Appendix F, gives a brief description of each subroutine's function and 
nesting. 

b. Math Models 

The math model for the UDU' factorization of the inverse of the 
B-matrix (which contains the partial derivatives of the observation equations with 
respect to the parameters being solved for) used in program KALMN2 is described in 
Bierman (1977) and Meyerho(T(1986). 

The TI-4100 receiver provides two pseudorange and two Doppler 
measurements (via LI and L2 frequencies), for as many as four satellites 
simultaneously. The following Equations (3.37 to 3.42) describe the various corrections 
and conversions; receiver clock bias, continuous count (Doppler) to interval (Doppler) 
data, and Ionospheric refraction, which are required for the Doppler and the 
pseudorange data (for the Ionospheric correction) (MeyerholT, 1985): 

Receiver Clock Bias Correction 

tt = ttr - RCB (3.37) 

where, 

ttr = receiver clocks GPS time-tag for range and delta 
range (Doppler) data 

tt = GPS time-tag corrected for receiver clock bias 

RCB = estimate of the receiver clock bias which is 

initially set to zero and an improved estimate 
of "RCB" is made after each time line 

The time delay data are multiplied by the speed of light "c" and 
converted to a pseudorange in [km]. 

Continuous Count to Interval ( Doppler ) Data Conversion 

DOP u (tt)j = DCj(tt)j - DCj(tt - td)j - (td * OFFSET1) (3.38) 
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DOP L2 (tt)j = DC 2 (tt) i - DC 2 (tt - td)j - (td * OFFSET2) 



(3.39) 



where, 

td = Doppler count interval for delta range (Doppler) data 



DCj(tt - td) , DC 2 (tt - td)= LI and L2 continuous 



Doppler count data at the delta range interval 

DCj(tt) , DC 2 (tt) = LI and L2 continuous count Doppler 
data at the end of the delta range interval 
OFFSET j , OFFSET 2 = Doppler offsets used by the receiver 
for LI and L2, -6000 and +7600 respectively 

DOPLj(tt) , DOP^ 2 (tt) = LI and L2 interval Doppler data 
for time-tag "tt" 

i = satellite from which data originates 

The delta range (Doppler) data is obtained by converting the LI interval Doppler data 
for each satellite. 



where, 

V os = nominal satellite oscillator frequency, 10.23 MHz 
Qj = LI frequency multiplier ( 154) 

RD(tt) = delta range (Doppler) for one satellite in [km] 
Ionospheric Refraction Correction 



RD(tt): = (-DOP ui (tt)/V 0S * Qj) * c 



(3.40) 



RCiaOj = CPI(tt)j + CIO(tt)i 



(3.41) 



RDCl(tt)j = RD(tt)j - DCIO(tt)i 



(3.42) 



where, 
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CPI = uncorrected pseudorange measurement in [km] 

RD(tt) = uncorrected delta (Doppler) range measurement 

CIO(tt) = ionospheric refraction correction in [km] 

for pseudorange data (refer to Equation 3.15) 

DCIO(tt) = ionospheric refraction corrections in [km] 

for delta (Doppler) .range data (refer to Equation 
3.16). 

Smoothing For Noise 

The pseudorangc data are contaminated due to signal multipath and 
other noise sources (MeyerhofT and Evans, 1986). To suppress the noise due to 
multipath in the pseudorange data, they are smoothed by the delta range (Doppler) 
data. The next set of Equations (3.43, 3.44, and 3.45) describes this smoothing: 

B N (tt) = R(tt) - DR(tt,ttfst) (3.43) 



B avg (tt) = lBj(tt)/N ' (3.44) 

R sm (tt) = B avg (tt) + DR(tt.ttfst) (3.45) 

where, 

tt = GPS time-tag for pseudorange and delta range (Doppler) data 

ttfst = GPS time-tag for the first good pseudorange and 
Doppler count of this data span 

R(tt) = Ionospheric corrected pseudorange in [km] 

DR(tt, ttfst) = Ionospheric corrected delta range (Doppler) data 

X = number of data points in the span 

13^- = estimate of pseudorange at time "ttfst" 
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B av „ = average value of the first pseudorange 

R sm (tt) = smooth pseudorange, if the absolute value of 

the difference between the raw pseudorange and 
smoothed pseudorange exceeds 20 [m], the 
data point is considered bad. 

Both the pseudorange and delta range (Doppler) data must be corrected for Time 

Epoch, clock errors, earth rotation, relativistic effects, and Tropospheric refraction. 

These formulae are the same ones used for the iterative least squares solution in the 

SOLTOM program (Equations 3.12, 3.13, 3.14, 3.31, 3.32, 3.33, 3.34, 3.35, 3.29, 3.30). 

The corrected pseudorange is given by, 



RC(tt)j = R sm (tt)i + C cl (tt)i - Rel(tt)j - [TR(tt) i * (1 + CR)J + C er (3.46) 

RDC(tt) is calculated in a similar manner as (3.46) for the delta range (Doppler) data. 

The following Equations (3.47 to 3.51) describe computed ranges and residuals for 
pseudorange and delta range (Doppler) data (Meyerhoff, 1985): 

Computed Ranges 

CR(tt)j = [(X(t)j - Xr) 2 + (Y(t)j-Yr) 2 + (Z(t)j - Zr ) 2 ] 1 / 2 (3.47) 

for the range at the beginning of the Doppler interval, 

CR(tt - td)j = [{X(t2) i - [Xr - (XVr * td)]} 2 + {Y(t2)j - [Yr - (YVr * td)]} 2 (3.48) 

+ {Z(t2)j - [Zr - (ZVr * td)[} 2 \ l l 2 

where, 

X(t), Y(t), Z(t) = position of i 1 * 1 satellite in an earth- 
fixed geocentric cartesian coordinate 
system at transit time t 

Xr, Yr, Zr = estimate of the receiver's position in an 

earth-fixed geocentric cartesian coordinate 
system 
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X(t2), Y(t2), Z(t2) = position of i^ satellite in an 

earth-fixed geocentric 
cartesian coordinate system at 
time t2 

XVr, YVr, ZVr = estimate of receiver's velocity 

tt = GPS time-tag corrected for receiver clock bias 

td = Doppler count interval for delta range (Doppler) 
data 

Therefore the computed range difference for the individual satellites is, 
CRD(tt); = CR(tt)j - CR(tt - td) ; 



Residuals 

The pseudorange residuals are given by, 

Romc(tt); = RC(tt); - CR(tt); - RCB 
and for the delta range (Doppler) data, 

RDomc(tt); = RDC(tt); - CRD(tt); - RCD 
where, 

RC(tt); = corrected pseudorange (Equation 3.46) 

RDC(tt); = corrected delta range (Doppler) 

CR(tt); = computed pseudorangc 

DCR(tt); = computed delta range (Doppler) 

RCB = estimate of the receiver clock bias 

RCD = estimate of frequency bias of the receiver oscillator 
which is also receiver clock drift. 



(3.49) 



(3.50) 



(3.51) 
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tt = GPS time-tag corrected for receiver clock bias 



td = Doppler count interval for delta range (Doppler) 
data 

i = i'-h satellite being tracked 

A state vector SP is corrected and improved over time, is composed of 
the receiver position in an earth-fixed geocentric cartesian coordinate system (Xr, Yr, 
Zr) and receiver clock bias (RCB). The last set of Equations (3.52 to 3.58) describes 
the observation equation and development of the state vector for the pseudorange data 
(Meyerhoff, 1985): 

Observation Equation 

D(tt)j = CR(tt)j + RCB (3.52) 

or rewriting it as, 

D(tt)j = [(X(t) i - Xr) 2 + (Y(t)j - Yr) 2 + (Z(t)j - Zr) 2 ] 1 / 2 + RCB (3.53) 

The H-matrix is defined as containing the partial derivatives of the 
observation equation with respect to the state vector: 

IIj i = dD(tt)/3Xr = -(X(t)j - Xr)/CR(tt)j (3.54) 

H 2 ’j = 0D(tt)/5Yr = -(Y(t)j - Yr)/CR(tt)j 

II 3 ’i = aD(tt)/aZr = -(Z(t)j - Zr)/CR(tt)j 

II 4)i = dD(tt)/5RCB = 1 

The PA-matrix is defined as a 4x4 variance-covariance matrix for 
position and receiver clock bias, with olT diagonal elements zeroed out and the diagonal 
elements set to the initial variance of Xr, Yr, Zr, and RCB. The Kalman vector K is 
given by the matrix equation, 

K(tt)j = PA(tt-) * I-l T j * [ Hj * PA(tt-) * H T j + Rj l' 1 (3.55) 



where, 
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PA(tt-) = variance-covariance matrix at time tt before it is 
updated for this data point 

R = variance on the pseudorange observation 

t Vi 

i = i in satellite being tracked 
Once the PA-matrix is updated it is defined by, 

PA(tt + ) = PA(tt-) - K(tt)j * Hj * PA(tt-) (3.56) 



State Vector 

The final two equations define the updated state vector SP(tt + ): 

SP(tt + ) = SP(tt-) + K(tt)j * Romc(tt)j (3.57) 

or 

Xr( + ) = Xr(-) + K^ttJj 1,1 Romc(tt) i (3.58) 

Yr( + ) = Yr(-) + K 2 (tt)j * Romc(tt) i 

Zr( + ) = Zr(-) + K^(tt)- * Romc(tt)j 

RCB( + ) = RCB(-) + K 4 ( tt)j * Romc(tt)j 

where ( + ) means the updated state vector and (-) means before the updated state 
vector. 

c. Input/ Output 

There are two input files for the KALMN2 program (Figure 3.1). The 
first file is an output file from the CON9TR program which contains records that are 
aligned with common time-tagged GPS data items. The second file is a user input file 
which sets the types of solutions and/or constraints on that solution, along with which 
data solutions will be printed out. The user input guide for this file is in (Appendix G). 

There are three possible output files for KALMN2 (Figure 3.1). The 
first file is the updated receiver position and time and any other print options that are 
set by the user input file. The second file is a position history file that is used in 
printing out the updated receiver position and time file. The third file is an optional 
GDOP file, that can be set in the user input file. 
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d. Program Modifications at NFS 

The source code from NSWC for the KALMN2 program was written in 
FORTRAN IV which has been modified to run at NPS on the IBM-3033 mainframe 
computer. Similar changes to the source code for the CON9TR and SOLTOM 
program were also made to the KALMN2 program including the executive files 
(Chapter III.A.2.). The original code was written to handle both the precise as well as 
the onboard broadcast ephemeris data in the navigation data message from the SV. 
The test data for KALMN2 was generated at NSWC using the broadcast ephemeris 
only, and as such, the modified program version at NPS has the precise ephemeris 
subroutines and corresponding, call, statements commented out. 

e. Solution Techniques Through Kalman Filter 

Using a UDL" factorization of an eight-state extended Kalman filter 
(with smoothed pseudorange measurements) in program KALMN2, there are nine 
types of solutions that may be calculated depending on the user input file (Appendix 
G). 

Note that, as was stated previously in Chapter III.B. 2. a., the velocity 
solutions from the four-velocity states are still being tested at NSWC. 

The initial receiver position (X, Y, Z) and clock biases are updated using 
the corrected state vector from the Kalman filter (Equations 3.57 and 3.58). If for a 
given time line the absolute value for the ionospheric correction CIO(tt) exceeds 20 m, 
that data is considered bad and Ragged as such, with a quality vector set to 15. Height 
constrained (especially in marine areas) allows the user's receiver to track four satellites 
with poor GDOP and still be able to produce a good solution. 
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IV. TEST RUN RESULTS AND VALIDATION 



The testing of the three programs CON9TR, SOLTOM, KALMN2 was 
accomplished by using two independent data sets in the SEF on 9-track magnetic tape 
as received from NSWC. Included also were hardcopy printouts, listing the solutions 
from those data sets using the CON9TR, SOLTOM, and KALMN2 programs on 
NSWC's CDC Cyber-865 computer. A direct comparison of the output listing files 
containing the receiver position solutions and times for the SOLTOM and KALMN2 
programs from both the NSWC's CDC Cyber-865 and the NTS IBM-3033 computers 
could then be made. 

After the original source codes and some algorithms were modified for 
compatibility on the NPS IBM-3033, the programs were test run, using the FORTVS 
compiler with an auto-double precision option discussed previously (Chapter III.A.2.). 

A. STATIC POSITION SOLUTION RESULTS 

Table II shows the differences between results computed at NTS and NSWC 
using the NSWC test data set (from Henderson Point, 85006) for the SOLTOM 
program. Running the CON9TR and SOLTOM programs, position solution results 
(X, Y, Z) were duplicated to the millimeter level for solution type seven; to the 
centimeter level for solution types one and two; to the decimeter level for solution 
types four and six; and to the meter level for solution type nine (Appendix E). 

The above differences in the SOLTOM solutions are considered insignificant due 
to the limitations of the observed data (Meyerhoff, 19S6). 

B. DYNAMIC POSITION SOLUTION RESULTS 

Table III shows the differences between results (beginning and ending solutions 
using the test data set) computed at NPS and NSWC using the NSWC test data (from 
Henderson Point, 95005) for the KALMN2 program. Running the CON9TR and 
KALMN2 programs, position solutions and clock biases were duplicated to the 
millimeter and 10'^ second level respectively. 

It should be noted that not all options (Chapter 1 1 LB. 2 and Appendix G) within 
the programs have been checked and tested. At this point it can be said that these 
programs have been established and validated on the IBM-3033 at NTS to the degree 
of being able duplicate NSWC's results if the same user options arc followed. 
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TABLE II 




CON9TR AND SOLTOM SOLUTION RESULTS 




Difference Between NPS and NSWC Solutions 


X,Y,Z 

Solution Position C 
Type (kilometers) ( 


,ock Bias 
seconds) 


Clock Drift 
( seconds) 


Number 

Of Data Points 


ONE 

Range Data 
Receiver Clock 


10~ 5 


10" 9 


10 -11 


1133 


TWO 

Doppler Data 
Receiver Clock 


10“ 5 


N/A 


10' 11 


1125 


FOUR 

Common Range 
Receiver Clock 


10“ 4 


N/A 


N/A 


832 


SIX 

Range Data 
SV Clock 


10“ 4 


10- 7 


10" 12 


1133 


SEVEN 

Doppler Data 
SV Clock 


10" 7 


N/A 


10" 11 


1125 


NINE 

Common Range 
SV Clock 


10~ 3 


10" 6 


10“ 12 


832 
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TABLE III 

CON9TR AND KALMN2 SOLUTION RESULTS 



Difference Between NPS and NSWC Solutions 



GPS Time-Tag 
(seconds) 


(X.Y,Z) 

(Kilometers) 


Clock Bias 
( seconds) 


Number of 
Doppler Edited 
Points 


First Solution 








14129.00 


10" 5 


10" 9 


N/A 


Last Solution 








17499. 00 


10" 5 


10" 6 


9 
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V. SUMMARY 



How to process the GPS satellite data from raw cassettes off the TI-4100 
GEOSTAR receiver to a reduced position and time solution has been discussed. 

The capability of processing raw data from cassettes to SEF at NPS is a task for 
the future. 

NPS at this time has the capability, starting with a SEF 9-track magnetic tape, to 
use the CON9TR program (where the common time-tagged data are aligned), to select 
either the SOLTO.Y1 (an iterative least squares for positioning a static receiver) or the 
KALMN2 program (which uses a Kalman filter for positioning a dynamic receiver) and 
to compute position and time solutions. The calculations for program CON9TR, and 
the math models for programs SOLTOM and KALMN2 were documented. The 
changes required to modify the original codes and algorithms for these programs 
written for NSWC's CDC Cyber-865 computer to make them compatible with the 
IBM-3033 computer at NPS have been described. Testing of the programs shows a 
successful conversion (with the restrictions discussed in Chapter IV.) of the programs 
to the IBM-3033 mainframe computer. This software forms a basis for a GPS 
programs library that will aid research using the GPS at NPS. 
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VI. RECOMMENDATIONS 



Since the conversion of the CON9TR, SOLTOM, and KALMN2 programs has 
been successfully tested using the two independent data sets from NSWC, it is 
recommended to explore other options within the programs set by the user input files 
for the SOLTOM and KALMN2 programs. Options that should be explored for the 
SOLTOM program include the use of the precise ephemeris and clock files, the use of 
more than one iteration for the solutions, and the calculation of Newton- Raphson 
range solutions followed by a comparison with the iterative batch least squares 
solutions for a given time line. The KALMN2 program should be explored by using 
the precise ephemeris and clock files, turning the height constraint on and off' in the 
user input file and comparing the solution results. The refraction correction should be 
turned on and off to see how the solutions are affected (this is another option in the 
user input file). 

It is further recommended to explore the possiblities of rewriting or adding codes 
and algorithms to improve the solution corrections for updating the four-state vectors 
of the velocity in the extended eight-state Kalman filter making them reliable for use in 
velocity solutions. 

In the future TI-4100 software for carrier phase data, differential positioning etc., 
should be acquired and modified to be compatible with the IBM-3033 computer at 
NPS. 
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APPENDIX A 

CON9TR SUBROUTINE DESCRIPTIONS AND NESTING 



INPUT - 


This routine reads the user inputs. 


GETDAT - 


This routine reads the next set of satellite 
tracking data from the 9-track Standard tape 
and also updates weather, Navigation data, 
almanac and calibration. 


RDDTL - 


This routine initializes the quality vector 
flags and controls branching of the program 
depending on a particular error flag. 


ER9TR - 


This routine controls branching of the .program 
for several error and terminating conditions. 


UDNTDT - 


This routine processes non-tracking data so that 
it can be used, after it determines which type 
of non-tracking data was read from the 9-track 
Standard tape. 


PRTSL - 


This routine prints tape solutions. 


UDWTDT - 


This routine gets the weather data from the 9-track 
Standard tape and stores it in common that will be 
used to calculate the tropospheric correction. 


UDSATD - 


This routine reads the 9 : track Standard tape and 
defines new satellite orbit clock parameters. 


FDTRN - 


This routine will give the tracker number for the 
satellite that the new data is for. 


DEDPOF - 


This routine stores the values for the LI and L2 
Doppler offsets so that they may be used in the 
program. 


PRTPDT - 


This routine prints any unused 9-track Standard 
tape data. 


UDTRD - 


This routine decides if the tracking data is 
continuous count Doppler range or data quality 
measurement blocks and then calls the correct 
subroutine to handle the correct data type. 


DEQUDT - 


This routine defines data qualitv and then stores 
the data quality information. The satellite 
tracker mode and the quality vectors arc used to 
determine if the data are good. The remaining data 
values in this block arc used to calculate errors 
on the tracking data. 


DERGDT - 


This routine stores the range data and the signal 
to noise data in correct data array that will Be 
processed. 


DEDPDT - 


This routine converts the accumulated Doppler 
phase data into range difference data. The 
Doppler phase dataNs corrected for offset, and 
then it. is differenced with preceding data if the 
preceding data was good. 
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GNINP - 


This routine reads the logical unit number of the 
input file and returns the information read one 
item at a time. 


GN1D - 


This. routine controls how data items are read bv 
reading a four digit data type number. 


1X9TR - 


This routine initializes the read packaee of the 
9-track Standard tape and sets up a data block 
for the data items read off of the tape. 


BLDATA - 


This routine contains the data that has been read 
off of the 9-track Standard tape. 


SDAT - 


This routine takes the site name and the data 
interval from the 9-track Standard tape read 
common and passes these values. 


FS1G - 


This routine determines the difference between 
the raw ranae and the range difference data for 
the satellite'"being tracked.^ 


BLOCK DATA ■ 


■ This routine was added at NFS to make this 
program compatible with the FORTRAN 77 
compilers on the IBM-3033 mainframe computer. 
It consists of common statements and data 
statements that the risorous compilers 
at NFS would not permit where tnev were placed 
in the original source code of CON9TR. 




MAIN PROGRAM COX9TR 
INPUT 
GETDAT 

RDDTL 

ER9TR 

UDNTDT 

PRTSL 

UDWTDT 

UDSATD 

FDTRN 

DEDPOF 

PRTPDT 

UDTRD 

DEQUDT 

FDTRN 

DERGDT 

FDTRN 

DEDPDT 

FDTRN 

GNINP 

GNID 

IN9TR 

BLDATA 

GNINP 

RDDTL 

SDAT 

FSIG 

BLOCK DATA 
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APPENDIX B 

INPUT GUIDE FOR PROGRAM CON9TR 



( by S. L. MeyerhofT, NSWC) 



This is the user input file for the program with variable names and descriptions. The 
input is format free. 

line l 

1: NSAT = number of trackers (4) 

2: TSS = GPS time of first data point to put on 

new file 



3: TES 

4: K2 
5: K4 
line 2 
1: X 

2: Y 

3: Z 



time to end data on new file, if EOF 
found before TES then data will end at 
EOF 

range variance bias factor (kilometer)^ 
doppler variance bias factor (kilometer)^ 



X coordinate of initial receiver position 
that will be put on the new files 

Y coordinate of initial receiver position 
that will be put on the new files 

Z coordinate of initial receiver position 
that will be put on the new files 
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APPENDIX C 

TEXAS INSTRUMENTS GPS SATELLITE RECEIVER FILE FORMAT 

( by S. L. Meyerhoff, 1984, XSWC ) 

Data from the DMA version of the TI-4100 GPS satellite receiver can be put 
into the following format. In this format all common time satellite tracking data will 
be in the same record block. The format for this file is as follows. 



var 


type 


meaning 


Record One 






word 1 ITYPE 


intg 


This is the data type for 






the next record 



Record Two 



If ITYPE equals one 


then record two will con 


tain initial information. 


word 1 RL10FF 


real 


frequency bias offset for LI 


word 2 RL20FF 


real 


frequency bias offset for L2 


word 3 PB 


real 


pressure in (millibars) 


word 4 TP 


real 


temperature (degrees Celsius) 


word 5 IID 


real 


relative humidity (percent) 


word 6 X 


real 


X coordinate of initial 






receiver position (kilometers) 


word 7 Y 


real 


Y coordinate of initial 






receiver position (kilometers) 


word 8 Z 


real 


Z coordinate of initial 






receiver position (kilometers) 


word 9 NTR 


intg 


maximum number of trackers used 






will be 4 or larger 


word 10 DTREC 


real 


user given data interval 






(seconds) 


word 11 SITE(l) 


char 


site identification for data 


word 12 SITE(2) 


char 


site identification for data 
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(continued) 



If ITYPE equals two then record two will contain satellite ephcmeris data (navigation 
data message) 



word 1 J 


intg 


tracker number 


word 2 ISAT(J) 


intg 


satellite prn number 


word 3 AS(1) 


real 


satellite clock bias in seconds 


word 4 AS(2) 


real 


satellite clock drift 






(seconds/ second) 


word 5 AS(3) 


real 


satellite clock aging 
(scconds/second*sccond) 


word 6 CICS 


real 


cross track cosine amplitude 
(radians) 


word 7 CISS 


real 


cross track sine amplitude 
(radians) 


word 8 CRCS 


real 


radial cosine amplitude 
(kilometers) 


word 9 CRSS 


real 


radial sine amplitude 
(kilometers) 


word 10 CUCS 


real 


in track cosine amplitude 
(radians) 


word 11 CUSS 


real 


in track sine amplitude 
(radians) 


word 12 DNS 


real 


corrections to mean motion 






(radians) 


word 13 CS 


real 


eccentricity 


word 14 I DOTS 


real 


inclination rate 






(radians/second) 


word 15 I OS 


real 


inclination (radians) 


word 16 MOS 


real 


mean anomaly (radians) 


word 17 OMEDS 


real 


right ascension rate 






(radians/second) 


word 18 OMEGS 


real 


right ascension (radians) 


word 19 SQAS 


real 


square root of semi-major 



axis (kilometers) 
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word 20 ADE 


real 


age of data cphemeris 
(seconds) 


word 21 TOCS 


real 


clock epoch (seconds) 


word 22 TOES 


real 


ephemeris epoch (seconds) 


word 23 IWKNOS 


intg 


week number 


word 24 WS 


real 


argument of perigee 
(radians) 


word 25 1EDATS 


intg 


sv health 


word 26 ADC 


real 


age of data clock (seconds) 



If ITYPE equals three then record two will contain the tracking information. 



word 1 TT 




real 


GPS time-tag for data 








that follows (seconds) 


word 2 ISAT(l) 




The PRN number for the 


to 


to 


intg 


satellite being tracked 


word NTR+ 1 ISAT(NTR) 




on trackers i to NTR 


word XTR + 2 


CR1(1) 




The LI pseudorange on 


to 


to 


real 


each tracker (seconds) 


word NTR*2 + 1 


CRl(NTR) 






word NTR*2 + 2 


CR2(1) 




The L2 pseudorange on 


to 


to 


real 


each tracker 


word NTR*3 + 1 


CR2(NTR) 




(seconds) 


word NTR*3 + 2 


DOPl(l) 




The LI continuous 


to 


to 


real 


count dopplcr data on 


word NTR*4 + 1 


DOPl(NTR) 




each tracker (cycles) 


word NTR*4 + 2 


DOP2( 1) 




The L2 continuous 


to 


to 


real 


count dopplcr data on 


word NTR*5 + 1 


DOP2(NTR) 




each tracker (cycles) 


word NTR*5 + 2 


SGR1(1) 




The variance on the LI 


to 


to 


real 


pseudorange 


word NTR*6+ 1 


SGRl(NTR) 




(kilometerkilo meter) 


word NTR*6 + 2 


SGR2(1) 




The variance on the L2 



53 



to 


to 


real 


pscudorangc 


word NTR*7+ 1 


SGR2(NTR) 




(kilometer*kilometer) 


word NTR*7 + 2 


SGDl(NTR) 




The variance on the LI 


to 


to 


real 


doppler 


word NTR*8 + 1 


SGDl(NTR) 




(kilometer*kilometer) 


word NTR*8 + 2 


SGD2( 1) 




The variance on the L2 


to 


to 


real 


doppler 


word NTR*9 + 1 


SGD2(NTR) 




(kilometer*kilometer) 


word NTR*9 + 2 


MQVEL(I) 




The quality vector for 


to 


to 


intg 


all trackers 


word NTR*10+ 1 


MQVEL(NTR) 




15 = loss of lock 
14 = bad data, no 
loss of lock 
1 = good range data 
bad doppler data 
0 = all data good 
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APPENDIX D 

SOLTOM SUBROUTINE DESCRIPTIONS AND NESTING 



CARDIO - 
LLXYZ - 

SIGXYZ - 

MATRIX - 
LIXV2F - 
VMULFF - 
VMULFM - 
INDAT - 
NEWSOL - 
LSTSOL - 

GETDAT - 
DEDPDT - 

DEFBR - 
CROAT - 

BSADR - 



CRDAT2 - 
ION - 
IONRD - 



This routine reads the input for program SOLTOM 
Irom the output of CON9TR. 

This routine converts position in the geodetic 
coordinate system to the earth-fixed geocentric 
coordinate system. 

This routine converts sigmas in the geodetic 
system into sigmas in the earth-fixe(T 
geocentric coordinate system. 

This routine performs matrix operations by 
calling routines in the IMSL library. 

This routine is in the IMSL and docs a matrix 
inversion. 

This routine is in the IMSL and does matrix 
multiplication. 

This routine is in the IMSL and does a 
transpose of matrix multiplication. 

This routine will pass initial values to some 
data arrays. 

This routine sets up the configuration for a 
for a given solution routine. 

This routine calculates a sequential least 
squares solution for the receiver position and 
clock. 

This routine reads the next set of satellite 
tracking data from CON9TR output file and it 
will also update weather, Navigation data, 
almanac and calibration. 

This routine converts the accumulated Doppler 
phase data into range difference data by 
correcting phase lor offset and then 
differencing it with data from before 
(if before data is good). 

This routine defines bias range 

This routine defines correction for receiver 
bias. 

This routine will correct data for receiver 
clock bias, drift and also correct GPS time 
tag for the end of week cross over. 

This routine converts data into range 
residuals and corrects data for error sources. 

This routine corrects for ionospheric 
correction. 

This routine corrects range difference data. 
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EPHDAT - 



CLOCK - 

SATPOS - 

ESTPOS - 
ZANG - 

RES - 

SMOT - 
TROP - 
REL - 
FSIG - 

SOLTNS - 
SMPT1 - 

UDRRC - 
SMPT2 - 

UDRDRC - 
SMPT3 - 



This routine defines the following values 
from the ephemeris data, satellite tracker 
number, satellite bias, drift and aging terms, 
semi-major axis of satellite orbit, 
eccentricity of the satellite orbit and epoch 
time for clock terms. 

This routine adds satellite and clock error 
and converts it to range measurements in 
kilometers and also performs the same for range 
difference data. 

This routine will find the position of the 
satellite from the Navigation Message or from 
the Precise Ephemeris. 

This routine finds the position of the 
satellite from the Precise Ephemeris data. 

This routine finds the Zenith angle by doting 
the X, Y, Z positions of the satellite receiver 
position vector to the vector from the receiver 
to the satellite. 

This routine finds the residual by calculating 
the range to a satellite and subtracting it 
from the observed range. This is 3150" 
performed for range difference data. 

This routine determines an earth rotation 
correction. 

This routine finds tropospheric refraction 
correction. 

This routine finds a relativistic correction 
to the data. 

This routine finds the variance (sigma squared) 
for the raw range and range difference data for 
the satellite being tracked. 

This routine calls the different routines to 
sum data in the correct part of the "B" matrix 
and "E" vector depending on which options have 
been set by the user for type of solution. 

This routine calls a routine to define additions 
to "B" matrix and "E" vector plus directing the 
called routine on where to store this 
information for the "B" matrix and "E" vector 
(for range data only). 

Update range solution matrix receiver clocks. 
Then store data in matrix so that solution can 
computed. 

This routine calls a routine to define 
additions to "B" matrix and "E" vector plus 
directing the called routine on where to store 
this information for "B" matrix and "E" vector 
(for range diifercnce data only) 

Update range diifercnce solution matrix 
receiver bias. This routine stores data in 
matrix so a solution can be computed. 

t his routine calls a routine to deline 
additions to "B" matrix and "E" vector, plus 
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NET RAP - 

CROUT - 
LINV3F - 
POSAVG - 
SMPT6 - 

UDRSC - 

BELIM - 
BMSTR - 
SMPT7 - 

UDRDSC - 
SMPT8 - 

SMPT9 - 

UDCRSC - 

ELIM - 
SMPT10 - 



directing the called routine on where to store 
this information for the "B" matrix and "E" 
vector (for range and ranee difference data 
only). 

This routine uses a Newton-Raphson method to 
calculate time bias adjustment and then improve 
the receiver initial position. 

This routine calls an IMSL routine that finds 
the determinate of a matrix. 

This is an IMSL routine that calculates the 
determinate of a matrix. 

This routine finds the average of all positions 
up to current time line. 

This routine calls a routine to define additions 
to "B" matrix and "E" vector plus directing the 
called routine on where to store this 
information (for range data only, satellite 
biases) 

This routine stores data in matrix so solution 
so solution can be computed after the ranee 
solution matrix and satellite clocks have been 
updated. 

This routine eliminates certain parameters in 
the "B" matrix. 

This routine stores data in a matrix so a 
solution can be computed. 

This routine calls a routine to define 
additions to. the "B" matrix and "E" vector 
plus directing the called routine on where to 
store this information (for range difference 
data only, satellite clock). 

This routine stores data in a matrix so that 
a solution can be computed after updating the 
range solution matrix and satellite bias. 

This routine calls a routine to define additions 
to "B" matrix and "E" vector plus directing the 
called routine where to store this information 
(for range and range difference data only). 

This routine calls a routine to define additions 
to the "B" matrix and "E" vector plus directing 
the called routine on where to store this 
information (for common time range data only). 

This routine stores data in a matrix so that 
a solution can be computed after updating; the 
common range solution matrix and satellite 
clocks. 

This routine finds data needed for elimination 
and then calls a routine to do that elimination. 

This routine calls a routine to define 
additions to the "B" matrix and "E" vector plus 
directing the called routine on where to store 
this information (for Doppler data used as bias 
range data). 
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BSRNG - 
NEWITR - 
BLOCK DATA 



This routine defines additions to the "B" 
matrix and "E" vector. 

This routine resets data so that a new iteration 
can be done. 

This routine was added at NPS to make this 
program compatible with the FORTRAN 77 
compilers on the IBM-3033 mainframe 
computer system. It contains common 
statements and data statements that 
NPS's compilers would not permit where 
they were placed in the original code of 
program SOLTOM. 



MAIN PROGRAM SOLTOM 
CARDIO 
LLLIXYZ 
SIGXYZ 
MATRIX 

LINV2F 
VMULFF 
VMULFM 
IN DAT 
NEWSOL 

LSTSOL 

MATRIX 

GETDAT 

DEDPDT 

DEFBR 

CRDAT 

BSADR 

CRDAT2 

ION 

IONRD 

EPIIDAT 

CLOCK 

SATPOS 

ESTPOS 

ZANG 

RES 

SMOT 

TROP 

REL 

FSIG 

SOLTNS 

SMPT1 

UDRRC 

SMPT2 

UDRDRC 

SMPT3 

UDRRC 
UDRDRC 
NET RAP 

CROUT 

LINV3F 

CRDAT 

POSAVG 

SMPT6 

UDRSC 

BELIM 

MATRIX 

LINV2F 

VMULFF 

VMULFM 

BMSTR 

SMPT7 
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UDRDSC 

BELIM 

BMSTR 

SMPT8 

UDRSC 

UDRDSC 

SMPT9 

UDCRSC 

ELIM 

BELIM 

BMSTR 

SMPT10 

BSRNG 

CRDAT 

UDRSC 



NEWSOL 

NEWITR 

GETDAT 
BLOCK DATA 
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APPENDIX E 

INPUT GUIDE AND SOLUTION TYPES FOR PROGRAM SOLTOM 



( by S. L. Mcyerhoff, NSWC) 



INPUT GUIDE 



Var 



Type Description 



CARD ONE Format (2F 10.0,615) 

(This line can oe read format free) 

col 1 to 10 

TSS real The GPS time tag for data at start 

of solution (seconds) 



col 1 1 to 20 
TES 



real 



The GPS time : tag of last data point 
to use in solutioif(seconds) 



col 21 to 25 (right justified) 

NPN integer The number of.time lines of data 

needed before first solution 



col 26 to 30 (right justified) 

NPB integer The number of time lines of data 

between solutions 



col 31 to 35 (right justified) 
ITRY integer 

col 36 to 40 (right justified) 
NSAT integer 

col 41 to 45 (right justified) 
IGT integer 



col 46 to 50 (right justified) 
ISP integer 



The maximum number of iterations 



The maximum number of trackers to 
process at once (range 1 to 4) 



= 1 if data source is Standard 9-track 
format 

= 2 if data source GESAR test file 
format 

= 3 if data source is NAV IF tape 



= 1 for satellite position from NAV 
data 

= 2 for satellite position from 
Celestial Earth fixed trajectory 
files 



col 51 to 55 (right justified) 

1AGE " integer = 1 solve for clock aging, otherwise 

not solved for 

col 56 to 60 (right justified) 

ICLOCK integer = 0 solve clock corrections arc from 

NAV data message 
= 1 clock corrections are post fit 
data 1'rom cards 



CARD TWO ( receiver position) Formal ( IIOJGIO.O ) 

(This card can be read format free) 



col 10 



60 



1FLAG 


integer 




= 0 for cartesian coordinates 
= 1 for geodetic coordinates 


col 11 to 20 
XL AT 


real 




latitude or X coordinate for position 
(degrees or kilometers) 


col 21 to 30 
YLOG 


real 




longitude or Y coordinate for position 
(degrees or kilometers) 


col 31 to 40 
ZHT 


real 




height or Z coordinate for position 
(kilometers) 


col 41 to 50 
TB 


real 




time bias for receiver clock 
(kilometers) 


col 51 to 60 
TD 


real 




time drift for receiver clock 
(kilometers/second) 


col 61 to 70 
TO 


real 




epoch for clock drift (given and 
solution) If negative, then first 
used time line is the epoch for the 
solution clock drift 


col 71 to 80 
AGE 


real 




aging for receiver clock 
(kilofnetcrs/sccond*second) 


CARD THREE (Reference Geoid) Format (2G 10.0) 

(This card can be read format free) 


col 1 to 10 
AG 


real 




semi-major axis of the reference geoid 
(kilometers) 


col 1 1 to 20 
ROBL 


real 




Obliquity of the reference geoid 


CARD FOUR ( Sigmas for Position) Format (7G10.0) 

(This card can be read format free) 


col 1 to 10 
SIGLT 


real 




sigma of X or latitude coordinate 
(kilometers) 


col 11 to 20 
S1GLG 


real 




sigma of Y or longitude coordinate 
(kilometers) 


col 21 to 30 
SIGHT 


real 




sigma ofZ or height coordinate 
(kilometers) 


col 31 to 40 
SIG(l) 


real 




sigma of refraction scaling constant 


col 41 to 50 
SIG(2) 


real 




sigma on receiver or satellite 1 time 
bias (kilometers) 


col 51 to 60 
SIG(3) 


real 




sigma on receiver or satellite 1 
clock drift (kilomctcrs/second) 
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col 61 to 70 

SIG(4) real sigma for aging term or satellite 2 

bias (kilometers/second*sccond) 

CARD FIVE ( Sigmas for Position ) Format (7G10.0) 

(This card can be read format free) 



col 1 to 10 
S1G(5) 


real 


siema on satellite 2 drift 
(kilo meters/ second) 


col 1 1 to 20 
S1G(6) 


real 


sigma for satellite 3 bias 
(kilometers) 


col 21 to 30 
S1G(7) 


real 


sigma for satellite 3 drift 
(kilometers/second) 


col 31 to 40 
SIG(8) 


real 


sigma for satellite 4 bias 
(kilometers) 


col 41 to 50 
S1G(9) 


real 


siema for satellite 4 drift 
(kilometers) 


SOLUTION TYPES 




CARD SIX ( Types of Solutions to do) Format (1615 ) 

(This card cannot be read format free) 


col 5 
NSOL(l) 


integer 


= 1 for range solution receiver clock 
0 = no 


col 10 
NSOL(2) 


integer 


= 1. for range difference solution 
receiver clock 
= 0 no 


col 15 
NSOL(3) 


integer 


= 1 for range and range difference 
solution receiver clock 
= 0 no 


col 20 
NSOL(4) 


integer 


= 1 for common time range solution 




receiver clock 


= 0 no 


col 25 
NSOL(5) 


integer 


= 1 for Ncwton-Raphson range solution 
= 0 no 


col 30 
NSOL(6) 


integer 


= 1 for range solution satellite 




clocks 


= 0 no 


col 35 
NSOL(7) 


integer 


= 1 for range difference solution 




satellite clock 


= 0 no 



col 40 
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NSOL(8) 


integer 


= 1 for range and range difference 




solution satellite clock 

= 0 no 


col 45 
NSOL(9) 


integer 


= 1 for common time range solution, 




satellite clock 


= 0 no 


col 50 
NSOL(IO) 


integer 


= 1 for doppler data used as bias 




range data 


= 0 no 


CARD SEVEN ( Variance Bias Factor ) Format ( 2G10.0 ) 

(This card can be read format free) 


col 1 to 10 
K2 


real 


bias factor for variance on range 
data 


col 1 1 to 20 
K4 


real 


bias factor for variance on doppler 
data 


CARD EIGHT ( Print Options ) Format (515) 

(This card can be read lormat free) 


col 5 
IPRT(l) 


integer 


= 1 print position solution (X,Y,Z) 
= 0 do not print 


col 10 
IPRT(2) 


integer 


= 1 print solution sigmas and clock 




terms 


= 0 do not print 


col 15 
1PRT(3) 


integer 


= 1 print data and time tags 
= 0 do not print 


col 20 
IPRT(4) 


integer 


= 1 print residuals, range, trop, 




data sigma 


= 0 do not print 


col 25 
1PRT(5) 


integer 


= 1 print satellite position 
= 0 do not print 


col 30 
IPRT(6) 


integer 


= 1 print residuals and all data 




corrections 


= 0 do not print 


col 35 
1PRT(7) 


integer 


= 1 print solution position in 




lat, long, height 


= 0 do not print 


col 40 
1PRT(8) 


integer 


= 1 print real time solutions from 
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9-track 



= 0 do not print 

CARD NINE ( Weather Data, Kiev. Cutoff) Format ( 4GI0.0 ) 

(This card can be read format free) 

col 1 to 10 



PB 


real 


pressure in millibars 


col 1 1 to 20 
TC 


real 


temperature in degrees Centigrade 


col 21 to 30 
HD 


real 


relative humidity in percent 


col 31 to 40 
EL 


real 


elevation angle cutoff (degrees) 


CARD TEN (Station Number ) Format (15) 

(This card can be read format free) 


col 1 to 5 
ISTA 


integer 


station identification 


CARD ELEVEN (Tracker I Clock Parameters ) Format (3G10.0) 


col 1 to 10 
AO 


real 


clock bias of satellite tracked on 
tracker 1 


col 11 to 20 

A1 


real 


clock drift of satellite tracked on 
tracker 1 


col 21 to 30 
A2 


real 


clock aging of satellite tracked on 
tracker! 



CARD TWELVE ( Tracker 2 Clock Parameters ) Format ( 3GI0.0 ) 

(This card can be read format free) 



col 1 to 10 



AO 


real 


clock bias of satellite tracked on 
tracker 2 


col 11 to 20 

A1 


real 


clock drift of satellite tracked on 
tracker 2 


col 21 to 30 
A2 


real 


clock aging of satellite tracked on 
tracker z 


CARD THIRTEEN (Tracker 3 Clock Parameters) Formal (3G10.0) 

(This card can be read format free) 


col 1 to 10 
AO 


real 


clock bias of satellite tracked on 
tracker 3 


col 11 to 20 

A1 


real 


clock drift of satellite tracked on 
tracker 3 


col 21 to 30 
A2 


real 


clock aging of satellite tracked on 
tracker 5 
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CARD FOURTEEN (Tracker 4 Clock Parameters) Format{3G 10.0) 

(This card can be read format free) 

col 1 to 10 



AO 


real 


clock bias of satellite tracked on 
tracker 4 


col 11 to 20 






A1 


real 


clock drift of satellite tracked on 
tracker 4 


col 21 to 30 






A2 


real 


clock aging of satellite tracked on 
tracker d 
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APPENDIX F 

KALMN2 SUBROUTINE DESCRIPTIONS AND NESTING 



CARDIO - 
LLHXYZ - 

SIGXYZ - 

MATRIX - 

LINV2F - 

VMULFF - 

VMULFM - 

INDAT - 

TIMEUP - 

PROMPT - 
GETDAT - 

NEWRAN - 

DEDPDT - 

DEFBR - 
CRDAT - 

BSADR - 

CRDAT2 - 

ION - 

IONRD - 
EPIIDAT - 

CLOCK - 



This routine reads input for program. 

This routine converts latitude, longitude and 
height to X, Y. Z in earth-frxed geocentric 
cartesian coordinates. 

This routine converts sigmas in the geodetic 
coordinate system into sigmas in the earth- 
fixed geocentric coordinate system. 

This routine performs matrix operations, by 
calling the proper routine in the IMSL library. 

This routine is in the IMSL and performs a 
matrix inversion. 

This routine is in the IMSL and performs a 
matrix multiplication. 

This. routine performs a transpose of 
multiplication of matrices. 

This routine gives the initial values to the 
"B" matrix (L in Kalman filter) 

This routine uses a UDU' factorization method 
to update state matrix. 

This routine aids in the UDU' factorization. 

This routine reads the next set of satellite 
tracking data from CON9TR output file. 

This routine reduces noise and multipath in the 
range data. 

This routine converts Doppler phase into range 
differenced data. 

This routine defines bias range. 

This routine defines corrections for receiver 
bias. 

This routine will correct data for receiver 
clock bias and drift and also correct GPS time- 
tag for the end of week cross over. 

This routine converts data into range residuals 
and corrects data for error sources. 

This routine corrects for ionospheric 
refraction. 

This routine corrects range difference data 

This routine defines the values from the 
ephemeris data. 

This routine adds satellite and clock error and 
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SATPOS - 


converts it to range measurements in kilometers 
and also perform the same for range difference 
data. 

This routine will find the. position of the 
satellite from the navigation data message or 
from.the Celestial Earth- Fixed Trajectory 
(Precise Ephemeris). 


ESTPOS - 


This routine finds the position of the 
satellite from the Precise Ephemeris data. 


ZANG - 


This routine finds the zenith angle by doting 
the X, Y, Z positions of the satellite 'receiver 
position vector to the vector from the receiver 
to the satellite. 


RES - 


This routine finds the residual for ranges. 


SMOT - 


This routine determines an earth rotation 
correction. 


TROP - 


This routine calculates a tropospheric 
refraction correction. 


REL - 


This routine finds a relativistic correction 
to the data. 


FSIG - 


This routine calculates the variance for raw 
pscudorange and range differenced data for the 
satellite being tracked. 


GDOP - 


This routine if turned on, writes a GDOP file. 


LINV3F - 


This routine is in the I.Y1SL and performs a 
matrix inversion 


KALFIL - 


This routine uses corrected data to do the 
Kalman filter. 


INKAL - 


This routine sets up initial covariance matrix 
in UDU' format. 


UDFACT - 


This routine puts Apriori covariance matrix 
into UDU' factorized matrix. 


REAL - 


This routine calculates partial derivatives of 
the data equation for pseudorange data. 


UPDATE - 


This routine will update the Kalman solution 
for one data point. 


DPKAL - 


This routine calculates partial derivatives of 
the data equation for delta (Doppler) range 
differences. 


CONST R - 


This routine will constrain. the solution to 
the surface of a sphere (height) and constrain 
velocity. 


XYZLLII - 


Converts X, Y, Z in an earth-fixed geocentric 
coordinate system to geodetic coordinate system. 


PRNSOL - 


This routine prints solutions. 


COVEVL - 


This routine finds the covariance matrix from 
the UDU' factorized matrix in U. 


BLOCK DATA - 


This routine was added at NPS to handle 
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data statements and common statements that 
had to be moved in the code or written, to 
be compatible with the IBM-3033 FORTRAN 77 
compilers. 



MAIN PROGRAM KALMN2 
CARDIO 

LLHXYZ 

S1GXYZ 

MATRIX 

LINV2F 

VMULFF 

VMULFM 

INDAT 

TIMEUP 

PROMPT 

GETDAT 

DEDPDT 

DEFBR 

CRD AT 

BSADR 

CRDAT2 

NEWRAN 

IONRD 

ION 
IONRD 
EP1IDAT 
CLOCK 2 
CLOCK 
SA I POS 

ESTPOS 

ZANG 

RES 

SMOT 

TROP 

REL 

FSIG 

GDOP 

LINV3F 

KALFIL 

INKAL 

UDFACT 

LLHXYZ 

RGKAL 

UPDATE 

DPKAL 

UPDATE 

CONSTR 

UPDATE 

XYZLLII 

PRNSOL 

COVEVL 

BLOCK DATA 
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APPENDIX G 

INPUT GUIDE FOR PROGRAM KALMN2 



Var Type 
CARD ONE 
col 1 to 10 

TSS real 



col 1 1 to 20 
TES real 



( by S. L. Meyerhof!*, NSWC) 

Description 

Format ( 2F10.0.6I5 ) 

(This card can be read format free) 

The GPS time-tag for data at start 
of solution (seconds) 



col 21 to 25 (right justified) 
NPN integer 

col 26 to 30 (right justified) 
NPB integer 

col 31 to 35 (right justified) 
ITRY *" integer 



col 36 to 40 (right justified) 
NS AT integer 



col 41 to 45 (right justified) 
IGT integer 



col 46 to 50 (right justified) 
ISP *" integer 



col 51 to 55 (right justified) 
ICLOCK integer 



col 56 to 60 (right justified) 
IAGE integer 



col 61 to 65 (right justified) 
IAGE integer 



The GPS time T tag of last data point 
to use in solution (seconds) 



The number of time lines between 
solutions 



The number of solutions between 
solutions being printed 

The number of solutions printed 
covariance matrix prints 

The number of trackers to process 
in the solution (range 1 to 4) 



= 1 if data source is standard 
9-track format 
= 2 if data source is other 



= 1 for satellite position from 
NAVIGATION DATA MESSAGE 
= 2 for satellite .position from 
Celestial Earth Fixed 
Trajectory files 



= 0 clock corrections are from 
NAVIGATION DATA MESSAGE 
= 1 clock corrections are post 
fit data from cards 



= 1 make GDOP file 
otherwise no file is made 



= 1 solve for clock aging 
otherwise no file is made 



CARD TWO {Receiver Position) Format (/ 10,7 G 10.0) 

( I his card can be read format free) 



col 10 
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I FLAG 


integer = 0 for Cartesian coordinates 

= I for Geodetic coordinates 


col 1 1 to 20 
XL AT 


real 


latitude or X coordinates for 
position (degrees or kilometers) 


col 21 to 30 
YLOG 


real 


longitude or Y coordinate for 
position (degrees or kilometers) 


col 31 to 40 
ZHT 


real 


height or Z coordinate for 
position (kilometers) 


col 41 to 50 
TB 


real 


time bias for receiver clock 
(kilometers) 


col 51 to 60 
TD 


real 


time drift for receiver clock 
(kilometers/second) 


col 61 to 70 
TO 


real 


epoch for clock drift (given and 
solution) if negative, then first 
used time line is the epoch for 
the solution clock drift 


col 71 to 80 
AGE 


real 


aging for receiver clock 
(Kilometers/second*second) 


CARD THREE ( Reference Geoid) Format ( 2GI0.0 ) 

(This card can be read format free) 

col I to 10 

AG real semi-major axis of the reference 

geoid (kilometers) 


col 1 1 to 20 
ROBL 


real 


obliquity of the reference geoid 


CARD FOUR ( Sigmas for 


Position) Formal (7 G 10.0) 

(This card can be read format free) 


col 1 to 10 
S1GLT 


real 


sigma of X or latitude coordinate 
(kilometer) 


col 11 to 20 
S1GLG 


real 


sigma of Y or longitude coordinate 
(kilometer) 


col 21 to 30 
SIGHT 


real 


sigma of Z or height coordinate 
(kilometers) 


col 31 to 40 
S1G(1) 


real 


sigma on refraction scaling constant 


col 41 to 50 
SIG(2) 


real 


sigma on receiver time bias 
(kilometer) 


col 51 to 60 
SIG(3) 


real 


sigma on receiver clock drift 
(kilometer/ second) 


col 61 to 70 
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S1G(4) real sigma of velocity in X direction 

(kilomcters/seco'nd) 

CARD FIVE {Sigmas for Velocity) format (7 G 10.0) 

(This card can be read format free) 



col 1 to 10 
SIG(5) 


real 


sigma of velocity in Y direction 
(kilomcters/second) 


col 11 to 20 
SIG(6) 


real 


sigma of velocity in Z direction 
(kilometers/second) 


col 21 to 30 
SIG(7) 


real 


sigma on receiver time bias 
(kilometers) 


col 31 to 40 
S1G(8) 


real 


sigma on receiver clock drift 
(kilomcters/second) 


col 41 to 50 
S1G(9) 


real 


sigma not used 


CARD SIX { Types of solutions to do) Format {1615) 

(This card cannot be read format free) 


col 5 
NSOL(l) 


integer 


= 1 for range position solution 
= 0 no 


col 10 
NSOL(2) 


integer 


= 1 Height constraint on position 
= 0 constraint only if less than 4 
satellites 


col 15 
NSOL(3) 


integer 


= 1 write position solution on file 
= 0 no solution saved 


col 20 
NSOL(4) 


integer 


= 1 to solve for refraction 
correction in position 
solution 
= 0 no 


col 25 
NSOL(5) 


integer 


= 1 to solve for refraction 
correction in velocity 
solution 
= 0 no 


col 30 
NSOL(6) 


integer 


= 1 to solve for clock drift 
from position solution 
= 0 no 


col 35 
NSOL(7) 


integer 


= 1 do a velocity solution 
(this has not been used 
before) 

= 0 no 


col 40 
NSOL(8) 


integer 


= 1 do a time update of 
covariance matrix 



= 0 no 
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integer 



col 45 



NSOL(9) 



= 1 not used 
= 0 no 



col 50 
NSOL(IO) 


integer 


= 1 for a constraint on 
velocity solution 






= 0 no 


CARD SEVEN ( Variance Bias Factor) Format ( 2G10.0 ) 

(This card can be read format free) 


col 1 to 10 
K2 


real 


bias factor for variance on range 
data 


col 1 1 to 20 
K4 


real 


bias factor for variance on 
Doppler data 



CARD EIGHT ( Print Options) Format (515) 

(This card cannot be read format free) 



IPRT(l) 


integer 


= 1 print position solution 
(XXZ) K . 

= 0 do not print 


col 10 
1PRT(2) 


integer 


= 1 print solution sigmas and 
clock terms 
= 0 do not print 


col 15 
IPRT(3) 


integer 


= 1 print data and time tags 
= 0 do not print 


col 20 
IPRT(4) 


integer 


= 1 residuals, range, trop, 
data sigma 
= 0 do not print 


col 25 
IPRT(5) 


integer 


= 1 print satellite position 
= 0 do not print 


col 30 
IPRT(6) 


integer 


= 1 residuals and all data 

corrections 

= 0 do not print 


col 35 
IPRT(7) 


integer 


= 1 print position in lat, long, 
and height 


col 40 
IPRT(8) 


integer 


= 1 print real time solutions from 
9-track 


CARD NINE ( 


weather I)ata s Flew Cutoff) Format ( 4GI0.0 ) 

( I his card can be read format free) 


col 1 to 10 
PB 


real 


pressure in millibars 


col 11 to 20 
TC 


real 


temperature in degrees Centigrade 
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col 21 to 30 



I ID 


real 




relative humidity in percent 


col 31 to 40 
EL 


real 




elevation angle cutoff (degrees) 


CARD TEN {Station Number) Format {Free) 

(This card can be read format free) 


col 1 to 10 

Qd) 


real 




time propagation noise for X 
velocity (kilometers) 


col 1 1 to 20 
Q(2) 


real 




time propagation noise for Y 
velocity (kilometers) 


col 21 to 30 
0(3) 


real 




time propagation noise for Z 
velocity (kilometers) 


col 31 to 40 

Q(4) 


real 




time propagation noise for time 
drift 


col 41 to 50 
S1GPS 


real 




sigma on position constraint 


col 51 to 60 
S1GVL 


real 




sigma on velocity constraint 


CARD ELEVEN {Diagonal elements of state transition matrix) 

(This card can be read format free) 

col 1 to 10 

PH 1( 1,1) real state transition matrix for 

parameter 1 


col 11 to 20 
PHI(2,2) 


real 




state transition matrix for 
parameter 2 


col 21 to 30 
PHI(3,3) 


real 




state transition matrix for 
parameter 3 


col 31 to 40 
PH1(4,4) 


real 




state transition matrix for 
parameter 4 


col 41 to 50 
PI 11(5,5) 


real 




state transition matrix for 
parameter 5 


col 51 to 60 
PH 1(6,6) 


real 




state transition matrix for 
parameter 6 


col 61 to 70 
PH 1(7,7) 


real 




state transition matrix for 
parameter 7 



CARD TWELVE ( Diagonal elements of slate transition matrix ) 

(This card can be read format free) 



col 1 to 10 

I’ll 1(8,8) real state transition matrix lor 

parameter 8 

col 1 1 to 20 
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PHI(9,9) 

CARD THIRTEEN 

col 1 to 10 
I IT 



real state transition matrix for 

parameter 9 

(Height for position constraint) 

( This card can be read format free) 



real Height above reference 

ellipsoid to constrain 
position solution to 
(kilometer) 
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APPENDIX H 

ABBREVIATIONS AND ACRONYMS 



GPS 


= 


Global Positioning System 


NPS 


= 


Naval Postgraduate School Monterey CA 


DMA 


= 


Defense Mapping Agency 


NOAA 


= 


National Oceanic and Atmospheric Administration 


NSVVC 


= 


Naval Surface Weapons Center 


sv 


= 


Space Vehicle 


SEP 


= 


Spherical Error Probability 


P-code 


= 


Precision code 


C/A-code 


= 


Coarse Acquisition code 


WGS 


= 


World Geodetic System 


PDOP 


= 


Position Dilution of Precision 


HDOP 


= 


Horizontal Dilution of Precision 


VDOP 


= 


Vertical Dilution of Precision 


TDOP 


= 


Time Dilution of Precision 


SEF 


= 


Standardized Exchange Format 


CIO 


= 


Conventional International Origin 


ARL 


= 


Applied Research Laboratories at the University ol 
at Austin 


IBM 


= 


International Business Machines 


CDC 


= 


Computer Data Corporation 


m 


= 


meters 


km 


= 


kilometers 


s 


= 


second 


Hz 


= 


Hertz 


MIIz 


= 


Megahertz 
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